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^^ I The implied volatility surface (IVS) is a fundamental building block in computational finance. We 

provide a survey of methodologies for constructing such surfaces. We also discuss various topics which 
influence the successful construction of IVS in practice: arbitrage-free conditions in both strike and time, 
Ph , how to perform extrapolation outside the core region, choice of calibrating functional and selection of 
s«»/ ' numerical optimization algorithms, volatility surface dynamics and asymptotics. 
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1 Introduction 

The geometric Brownian motion dynamics used by Black and Scholes (1973) and Merton (1973) to price 
options constitutes a landmark in the development of modern quantitative finance. Although it is widely 
acknowledged that the assumptions underlying the Black-Scholes-Merton model (denoted BSM for the 
rest of the paper) are far from realistic, the BSAi formula remains popular with practitioners, for whom 
it serves as a convenient mapping device from the space of option prices to a single real number called 
the implied volatility (IV). This mapping of prices to implied volatilities allows for easier comparison of 
options prices across various strikes, expiries and underlying assets. 

When the implied volatilities are plotted against the strike price at a fixed maturity, one often observes a 
skew or smile pattern, which has been shown to be directly related to the conditional non-normality of the 
underlying return risk-neutral distribution. In particular, a smile reflects fat tails in the return distribution 
whereas a skew indicates return distribution asymmetry. Furthermore, how the implied volatility smile 
varies across option maturity and calendar time reveals how the conditional return distribution non- 
normality varies across different conditioning horizons and over different time periods. For a fixed relative 
strike across several expiries one speaks of the term structure of the implied volatility. 

We mention a few complications which arise in the presence of smile. Arbitrage may exist among the 
quoted options. Even if the original market data set does not have arbitrage, the constructed volatility 
surface may not be arbitrage free. The trading desks need to price European options for strikes and 
maturities not quoted in the market, as well as pricing and hedging more exotic options by taking the 
smile into account. 

Therefore there are several practical reasons [G2] to have a smooth and well-behaved implied volatility 
surface (IVS): 

1. market makers quote options for strike-expiry pairs which are illiquid or not listed; 

2. pricing engines, which are used to price exotic options and which are based on far more realistic 
assumptions than BSA4 model, are calibrated against an observed IVS; 

3. the IVS given by a listed market serves as the market of primary hedging instruments against 
volatility and gamma risk (second-order sensitivity with respect to the spot); 

4. risk managers use stress scenarios defined on the IVS to visualize and quantify the risk inherent to 
option portfolios. 

The IVS is constructed using a discrete set of market data (implied volatilities or prices) for different 
strikes and maturities. Typical approaches used by financial institutions are based on: 

• (local) stochastic volatility models 

• Levy processes (including jump diffusion models) 

• direct modeling of dynamics of the implied volatility 

• parametric or semi-parametric representations 

• specialized interpolation methodologies 

Arbitrage conditions may be implicitly or explicitly embedded in the procedure 

This paper gives an overview of such approaches, describes characteristics of volatility surfaces and 
provides practical details for construction of IVS. 



2 Volatility surfaces based on (local) stochastic volatility models 

A widely used methodology employs formulae based from stochastic volatility models to fit the set of 
given market data. The result is an arbitrage free procedure to interpolate the implied volatility surface. 
The most commonly considered stochastic volatility models are Heston and SABR and their extensions 
(such as time dependent parameters, etc) and we will concentrate on these models as well. Having time 
dependent parameters allows us to perform calibration in both strike and time directions. This is arguably 
better than the case of using constant parameter models in capturing inter-dependencies of different time 
periods. The main disadvantage when using time dependent parameters is the increase in computational 
time, since in many cases we do not have analytical solutions/approximations and we have to resort to 
numerical solutions when performing the calibration. However, for the considered time dependent models, 
namely Heston and SABR, (semi) analytical approximations are available, which mitigates this issue. 

We will also consider the hybrid local stochastic volatility models, which are increasingly being preferred 
by practitioners, and describe efficient calibration procedures for such models. 

2.1 Heston model and its extensions 

The Heston model is a lognormal model where the square of volatility follows a Cox-Ingersoll-Ross (CIR) 
process. The call (and put) price has a closed formula through to a Fourier inversion of the characteristic 
function. Details on efficient ways to compute those formulas are given in [98], recent advances were 
presented in [9-3], while [79] contains details for improving the numerical calibration. 

It is our experience, confirmed by discussions with other practitioners, that having only one set of pa- 
rameters is usually not enough to match well market data corresponding to the whole range of expiries used 
for calibration, especially for FX markets. Consequently we need to consider time dependent parameters. 

When the parameters are piecewise constant, one can still derive a recursive closed formula using a 
PDE/ODE methodology [114] or a Markov argument in combination with affine models [59], but formula 
evaluation becomes increasingly time consuming. 

A better and more general approach is presented in [IS], which is based on expanding the price with 
respect to the volatility of volatility (which is quite small in practice) and then computing the correction 
terms using Malliavin calculus. The resulting formula is a sum of two terms: the BSM price for the model 
without volatility of volatility, and a correction term that is a combination of Greeks of the leading term 
with explicit weights depending only on the model parameters. 

The model is defined as 



with initial conditions 
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Using the same notations as in [IS], the price for the put option is 



Put{K, T) = exp 
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where r(t),q{t) are the risk free rate and, respectively, dividend yield, K is the strike and T is the 
maturity of the option. 

There are two assumptions employed in the paper: 

1) Parameters of the CIR process verify the following conditions 
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2) Correlation is bounded away from -1 and 1 
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Under these assumptions, the formula for approximation is 
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where PBs{x,y) is the price in a BSM model with spot e^, strike K, total variance y, maturity T and 
rates req,qeq given by 
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while var{T),ai{T),b2i{T) have the following formulas 
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The error is shown to be of order O {Csup^^) 

While results are presented for general case, [18] also includes explicit formulas for the case of piecewise 
constant parameters. Numerical results show that calibration is quite effective and that the approximation 
matches well the analytical solution, which requires numerical integration. They report that the use of 
the approximation formula enables a speed up of the valuation (and thus the calibration) by a factor 100 
to 600. 

We should also mention the efficient numerical approach presented in [] JO] for calibration of the time 
dependent Heston model. The constrained optimization problem is solved with an optimization procedure 
that combines Gauss-Newton approximation of the Hessian with a feasible point trust region SQP (sequen- 
tial quadratic programming) technique developed in [146]. As discussed in a later chapter on numerical 
remarks for calibration, in the case of local minimizer applied to problems with multiple local/global 
minima, a regularization term has to be added in order to ensure robustness/stability of the solution. 

2.2 SABR model and its extensions 

The SABR model [So] assumes that the forward asset price F{t) and its instantaneous volatility a{t) are 
driven by the following system of SDEs: 

dF{t) = a{t)F'^{t)dWi{t) (2.3) 

da{t) = iya{t)dW2(t) 
d{Wi,W2) = pdt 

where is z^ > is volatility of volatility and /3 > is a leverage coefficient. The initial conditions are 

F{0) = Fo 
a(0) = ao 

Financial interpretation for this model is the following: a{t) determines the overall level of at-the-money 
forward volatility; /? measures skew with two particular choices: /3 = 1 corresponding to the log-normal 
model with a flat skew and /3 = corresponding to the normal model with a very pronounced skew; p 
also controls the skew shape with the choice p < (respectively p > 0) yielding the negative (respectively 
inverse) skew and with the choice p = producing a symmetric volatility smile given /3 = 1; i/ is a measure 
of convexity, i.e. stochasticity of a{t). 

Essentially, this model assumes CEV distribution (log-normal in case /3 = 1) for forward price F{t) and 
log-normal distribution for instantaneous volatility a{t). 

SABR model is widely used by practitioners, due to the fact that it has available analytical approxi- 
mations. Several approaches were used in the literature for obtaining such approximations: the singular 
perturbation, heat kernel asymptotics, and Malliavin calculus [94, 115, 85]. Additional higher order ap- 
proximations are discussed in [119](second order) and [134], up to fifth order. Details for improving the 
numerical calibration were given in [79] 

An extension of SABR (termed lambda-SABR), and corresponding asymptotic approximations were 
introduced in papers by Henry-Labordere (see chapter 6 of [87]). This model is described as follows (and 
degenerates into SABR model when A = 0) 

dF{t) = a{t)F^{t)dWi{t) (2.4) 

da{t) = A (a(t) - X) + ua{t)dW2{t) 
d{Wi,W2) = pdt 



The high order approximations in [134] were obtained for lambda-SABR model. Approximations for 
extended lambda-SABR model, where a drift term is added to the SDE describing F{t) in (2.4), were 
presented in [130] and [131]. Approximations for SABR with time dependent coefficients were presented 
in [118], where the model was named "Dynamic SABR", and in respectively in [Nd], where the approach 
was specialized to piecewise constant parameters. 

If one combines the results from [131, 134, 130] with the findings of [80], the result will be a model 
(extended lambda-SABR with piecewise constant parameters) that may be rich to capture all desired 
properties when constructing a volatility surface and yet tractable enough due to analytical approxima- 
tions. 

Alternatively, the results presented in [80] seem very promising and will be briefly described below. The 
procedure is based on asymptotic expansion of the bivariate transition density [147]. 

To simplify the notations, the set of SABR parameters is denoted hy 9 = {a, /3, p, i'} and the dependence 

of the model's joint transition density on the model parameters hy p (t, Fq, aoj T, F,A; 6 1 . 
The joint transition density is defined as 

f(f < F{T) <F + dF,A< a{T) <A + dA\ =p ft, F, A; T, F, A) dPdA 

We follow the notations from [i 17], namely: 

• F,A are forward variables denoting the state values of F{T),a{T) 

• F,A are backward variables denoting the state values of F{t), a{t) 

Let us denote by {Ti, r2, ..., T/v} the set of expiries for which we have market data we want to calibrate to; 
we assume that the four SABR parameters {a,f3,p.i'} are piecewise constant on each interval [T!j_i,Tj]. 
The tenor-dependent SABR model then reads 

(2.5) 
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(2.5) is considered together with 
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The notations for SABR set of parameters and, respectively, for the dependence of the model's joint 
transition density on the model parameters are updated as follows 



9{T) 
and, respectively 
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(ai_i,/3j_i,pj_i,fj_i) ifTi_i <T <Ti 
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In the case where only the parameter dependence need to be stressed, we use the shortened notion: 



p(O,ri;0o) 
p(ri_i,Ti;6'i_i) 

A standard SABR model describes the dynamics of a forward price process F (t;Ti) maturing at a 
particular Tj. Forward prices associated with different maturities are martingales with respect to different 
forward measures defined by different zero-coupon bonds B {t,Ti) as numeraires. This raises consistency 
issues, on both the underlying and the pricing measure, when we work with multiple option maturities 
simultaneously. 

We address this issue by consolidating all dynamics into those of F (t, T^) , a (t, T/v), whose tenor is the 
longest among all, and express all option prices at different tenors in one terminal measure Q'^^ which is 
the one associated with the zero-coupon bond B (t,T/v) • 

We may do so because we assume 

• No-arbitrage between spot price S{t) and all of its forward prices F {t,Ti) ,i = 1...N, at all trading 
time t; 

• Zero-coupon bonds B {t,Ti) are risk-less assets with positive values 
Based on these assumptions we obtain the following formulas 
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This will enable us to convert an option on F{-,Ti) into an option on F{-,Tis[). The price of a call 
option on F {-^Ti) with strike price Kj and maturity Tj then becomes 



where the adjusted strike Kj is defined as 
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In the context of model calibration, computation of spot implied volatilities from the model relies on 
computation of option prices 



[F{T,)-Kjy\F,^uAi-i 
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dFidAi 



(2.6) 



at each tenor Tj for each equivalent strike Kj . 

Asymptotic expansions of a more generic joint transition density have been obtained analytically in 
[147] to the n*^ order. We should note that for simplicity we use exactly the same notations as in section 
4.2 of [80], namely that the values of state variables at time t are denoted by /, q and, respectively, the 
values of the state variables at T are denoted by F,A. 

The expansion to second order was shown to give a quite accurate approximation 
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The terms Po,pi,p2 have the fohowing formulae 
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Exphcit expressions for the polynomial functions an, aio, a23) 022j 02I) Q^20 ^-re given in Eq. (42) in [147]. 
In terms of computational cost, it is reported in [SO] that it takes about 1-10 milliseconds for an evaluation 
of the integral (2.6) on a 1000 by 1000 grid, using the approximation p2 as density. 

2.3 Local stochastic volatility model 

More and more practitioners are combining the strengths of local and stochastic volatility models, with 
the resulting hybrid termed local stochastic volatility (LSV) model. 

Based on [J 10], we describe efficient procedures for calibrating one such model, namely a hybrid Heston 
plus local volatility model, with dynamics given by 



a(/^^^(t),t)v^/^^^(t)dVFi(t) 
k{9- v{t)) dt + cVW)dW2{t) 



dv{t) 
The calibration procedure is based on the following 2 step process: 

• calibrate stochastic volatility component 

• perform LSV correction 

The validity of this 2-step process is due to the observation that the forward skew dynamics in stochastic 
volatility setting are mainly preserved under the LSV correction. 

The first approach is based on the "fixed point" concept described in [126] 

1. Solve forward Kolmogorov PDE (in x = In (•S'//wd) with a given estimate of cr{f, t) 



dp 
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2. Use the density from 1. to compute the conditional expected value of v{t) given f^^^{t) 

Jo P{t,f,v)dv 

3. Adjust a according to Gyongy's identity [S3] for the local volatilities of the LSV model 

KD' if,t) = <TM,t)E [v{t)\f'^{t) = f] = {atr'^^f {f,t) 

4. repeat steps 1.-3. until (T{f,t) has converged (it was reported that in most cases 1-2 loops are 
sufficient) 

The second approach is based on "local volatility" ratios, similar to [120, 88]. The main idea is the following: 
applying Gyongy's theorem [S-]] twice (for the starting stochastic volatility component and, respectively, 
for the target LSV model) avoids the need for conditional expectations. 
The procedure is as follows 

1. Compute the local volatilities of an LSV and an SV model via Gyongy's formula 

(^£r)' (/,*) = '^'(ME Hm^'^'it) = f] = (affy-'^f {f,t) 
{al^)\x,t) = E[v{t)\f''{t) = x] 

2. Taking the ratio and solving for the unknown function a{-, ■) we obtain 



"^''^^ alU-,t) V^Kt)l/^^^(t) = /] "''"^" ""^^''^ crlUHif,t),t) 

with an approximate, strictly monotonically increasing map H{f, t) 

The calculation is reported to be extremely fast if the starting local volatilities are easy to compute. The 
resulting calibration leads to near perfect fit of the market 

We should also mention a different calibration procedure for a hybrid Heston plus local volatility model, 
presented in [60]. 

3 Volatility surfaces based on Levy processes 

Volatility surface representations based on Levy processes tend to better handle steep short term skews 
(observed especially in FX and commodity markets). In a model with continuous paths like a diffusion 
model, the price process behaves locally like a Brownian motion and the probability that the price of 
the underlying moves by a large amount over a short period of time is very small, unless one fixes an 
unrealistically high value of volatility. Thus in such models the prices of short term OTM options are 
much lower than what one observes in real markets. On the other hand, if price of underlying is allowed 
to jump, even when the time to maturity is very short, there is a non-negligible probability that after a 
sudden change in the price of the underlying the option will move in the money. 
The Levy processes can be broadly divided into 2 main categories: 
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• jump diffusion processes: jumps are considered rare events, and in any given finite interval there are 
only finite many jumps 

• infinite activity Levy processes: in any finite time interval there are infinitely many jumps. 

The importance of a jump component when pricing options close to maturity is also pointed out in the 
literature, e.g., [4]. Using implied volatility surface asymptotics, the results from [111] confirm the presence 
of a jump component when looking at S&P option data . 

Before choosing a particular parametrization, one must determine the qualitative features of the model. 
In the context of Levy-based models, the important questions are [42, 138]: 

• Is the model a pure-jump process, a pure diffusion process or a combination of both? 

• Is the jump part a compound Poisson process or, respectively, an infinite intensity Levy process? 

• Is the data adequately described by a time-homogeneous Levy process or is a more general model 
may be required? 

Well known models based on Levy processes include Variance Gamma [JOT], Kou [100], Normal Inverse 
Gaussian [12], Meixner [129, 108], GGMY [33], affine jump diffusions [55]. 

From a practical point of view, calibration of Levy-based models is definitely more challenging, especially 
since it was shown in [43, 44] that it is not sufficient to consider only time- homogeneous Levy specifications. 
Using a non-parametric calibration procedure, these papers have shown that Levy processes reproduce the 
implied volatility smile for a single maturity quite well, but when it comes to calibrating several maturities 
at the same time, the calibration by Levy processes becomes much less precise. Thus successful calibration 
procedures would have to be based on models with more complex characteristics. 

To calibrate a jump-diffusion model to options of several maturities at the same time, the model must 
have a sufficient number of degrees of freedom to reproduce different term structures. This was demon- 
strated in [139] using the Bates model, for which the smile for short maturities is explained by the presence 
of jumps whereas the smile for longer maturities and the term structure of implied volatility is taken into 
account using the stochastic volatility process. 

In [74] a stochastic volatility jump diffusion model is calibrated to the whole market implied volatility 
surface at any given time, relying on the asymptotic behavior of the moments of the underlying distribution. 
A forward FIDE (Partial Integro-Differential Equation) for the evolution of call option prices as functions 
of strike and maturity was used in [4] in an efficient calibration to market quoted option prices, in the 
context of adding Poisson jumps to a local volatility model. 

3.1 Implied Levy volatility 

An interesting concept was introduced in [45], which introduced the implied Levy space volatility and the 
implied Levy time volatility, and showed that both implied Levy volatilities would allow an exact fit of 
the market. Instead of normal distribution, as is the case for implied volatility calculation, their starting 
point is a distribution that was more in line with the empirical observations. 
More specifically, instead of lognormal model they assume the following model 

S{t) = Soexp[{r-q + uj)t + aX{t)] (3.1) 

where cr > 0, r is the risk-free rate, q is the dividend yield, a; is a term that is added in order to make 
dynamics risk neutral, and X = {X{t),t > 0} is a stochastic process that starts at zero, has stationary 
and independent increments distributed according to the newly selected distribution 
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Once one has fixed the distribution of X (which we assume as in the Brownian case to have mean zero 
and variance at unit time equal to 1), for a given market price one can look for the corresponding o", which 
is termed the implied (space) Levy volatility, such that the model price matches exactly the market price. 

To define Implied Levy Space Volatility, we start with an infinitely divisible distribution with zero mean 
and variance one and denote the corresponding Levy process by X = {X{t),t > 0} . Hence £J[X(1)] = 
and yar[X(l)]=0. We denote the characteristic function of X{1) (the mother distribution) by 4'{u) = 
E[exp{iuX {!))]. We note that, similar to a Brownian Motion, we have £[X(i)] = 1 and yar[X(i)]=t and 
hence Var[aX{t)] = aH. 

If we set w in (3.1) to be a; = — log ((/>(— cii)), we call the volatility parameter a needed to match the 
model price with a given market price the implied Levy space volatility of the option. 

To define the implied Levy time volatility we start from a similar Levy process X and we consider that 
the dynamics of the underlying are given as 

S{t) = 5oexp[(r-g + wa2)t + X(i)] (3.2) 

cj = log (<?;>(-i)) 

Given a market price, we now use the terminology of implied Levy time volatility of the option to 
describe the volatility parameter a needed to match the model price with the market price. Note that 
in the BSM setting the distribution (and hence also the corresponding vanilla option prices) of aW{t) 
and W{a'^t) is the same, namely a A/'(0, a'^t) distribution, but this is not necessary the case for the more 
general Levy cases. 

The price of an European option is done using characteristic functions through the Carr-Madan formula 
[35] and the procedure is specialized to various Levy processes: normal inverse Gaussian (NIG), Meixner, 
etc. 

4 Volatility surface based on models for the dynamics of implied 
volatility 

In the literature there are two distinct directions for treatment and construction of volatility surfaces 
[36]. One approach assumes dynamics for the underlying that can accommodate the observed implied 
volatility smiles and skews. As we have seen in previous chapters, such approaches include stochastic 
volatility models as well as various Levy processes. The general procedure is to estimate the coefficients 
of the dynamics through calibration from observed option prices. Another approach relies on explicitly 
specifying dynamics of the implied volatilities, with models belonging to this class being termed "market 
models" of implied volatility. In general, this approach assumes that the entire implied volatility surface 
has known initial value and evolves continuously over time. The approach first specifies the continuous 
martingale component of the volatility surface, and then derives the restriction on the risk-neutral drift 
of the surface imposed by the requirement that all discounted asset prices be martingales. Such models 
are presented in [9, 61, 84, 48] An approach that was described as falling between the two categories was 
described in [361 and is described next 



4.1 Carr and Wu approach 

Similar to the market model approach, it directly models the dynamics of implied volatilities. However, 
instead of taking the initial implied volatility surface as given and infer the risk-neutral drift of the implied 
volatility dynamics, both the risk-neutral drift and the martingale component of the implied volatility 
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dynamics are specified, from which the ahowable shape for the implied volatility surface is derived. The 
shape of the initial implied volatility surface is guaranteed to be consistent with the specified implied 
volatility dynamics and, in this sense, this approach is similar to the first category. 

The starting point is the assumption that a single standard Brownian motion drives the whole volatility 
surface, and that a second partially correlated standard Brownian motion drives the underlying security 
price dynamics. By enforcing the condition that the discounted prices of options and their underlying are 
martingales under the risk-neutral measure, one obtains a partial differential equation (PDE) that governs 
the shape of the implied volatility surface, termed as Vega-Gamma- Vanna- Volga (VGVV) methodology, 
since it links the theta of the options and their four Greeks. Plugging in the analytical solutions for the 
BSM Greeks, the PDE is reduced into an algebraic relation that links the shape of the implied volatility 
surface to its risk- neutral dynamics. 

By parameterizing the implied variance dynamics as a mean-reverting square-root process, the algebraic 
equation simplifies into a quadratic equation of the implied volatility surface as a function of a standardized 
moneyness measure and time to maturity. The coefficients of the quadratic equation are governed by six 
coefficients related to the dynamics of the stock price and implied variance. This model is denoted as the 
square root variance (SRV) model. 

Alternatively, if the implied variance dynamics is parametrized as a mean-reverting lognormal process, 
one obtains another quadratic equation of the implied variance as a function of log strike over spot and 
time to maturity. The whole implied variance surface is again determined by six coefficients related to the 
stock price and implied variance dynamics. This model is labeled as the lognormal variance (LNV) model. 

The computational cost for calibration is quite small, since computing implied volatilities from each of 
the two models (SRV and LNV) is essentially based on solving a corresponding quadratic equation. 

The calibration is based on setting up a state-space framework by considering the model coefficients 
as hidden states and regarding the finite number of implied volatility observations as noisy observations. 
The coefficients are inferred from the observations using an unscented Kalman filter. 

Let us introduce the framework now. We note that zero rates are assumed without loss of generality. 

The dynamics of the stock price of the underlying are assumed to be given by 



dS{t) = S{t)y^dW{t) 

with dynamics of the instantaneous return variance v{t) left unspecified. 

For each option struck at K and expiring at T, its implied volatility I{t; K, T) follows a continuous 
process given by 

dl{t] K, T) = ^{t)dt + u{t)dZ{t) 

where Z{t) is a Brownian motion. 

The drift fi{t) and volatility of volatility uj{t) can depend on K, T and I(t; K, T) 

We also assume that we have the following correlation relationship 

p{t)dt = E [dW{t)dZ{t)] 

The relationship /(t; K,T) > guarantees that there is no static arbitrage between any option at {K; T) 
and the underlying stock and cash. 

It is further required that no dynamic arbitrage (NDA) be allowed between any option at (K; T) and 
respectively a basis option at [Kq]Tq) and the stock. 

For concreteness, let the basis option be a call with C{t\ T, K) denoting its value, and let all other 
options be puts, with P(t; K, T) denoting the corresponding values. We can write both the basis call and 
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other put options in terms of the BSAi put formula: 

P{t;K,T) = BSM{S{t),I{t;K,T),t) 
C{t;Ko,To) = BSM{S{t),I{t;Ko,To),t) + S{t)-Ko 

We can form a portfoho between the two to neutrahze the exposure on the volatihty risk dZ 
^^^ (S{t),I {t;K,T) ,t)u{K,T) - N'{t)^^^ {S{t),I {t;Ko,To) ,t)uj{Ko,To) = 

We can further use N (t) shares of the underlying stock to achieve delta neutrality: 

BSM {S{t),I (t; K, T),t)- iV^(t) [1 + BSM {S{t),I (t; Kq^Tq) , t)] - N^{t) = 

Since shares have no Vega, this three-asset portfolio retains zero exposure to dZ and by construction 
has zero exposure to dW. 

By Ito's lemma, each option in this portfolio has risk-neutral drift (RND) given by 

„,,^ dBSM , ,dBSM v(t) ^o, .d'^BSM ,, ,, r^^^, .d'^BSM co'^ (t) d'^ B S M 

Note: For simplicity of notations, for the remainder of the chapter we will use B instead of BSM 
No dynamic arbitrage and no rates imply that both option drifts must vanish, leading to the fundamental 
"PDE". 

The class of implied volatility surfaces defined by the fundamental "PDE" (4.1) is termed the Vega- 
Gamma- Vanna- Volga (VGVV) model 

We should note that (4.1) is not a PDE in the traditional sense for the following reasons 

• Traditionally, a PDE is specified to solve the value function. In our case, the value function 
B {S{t),I{t; K, T), i) is well-known. 

• The coefficients are deterministic in traditional PDEs, but are stochastic in (4.1) 

The "PDE" is not derived to solve the value function, but rather it is used to show that the various 
stochastic quantities have to satisfy this particular relation to exclude dynamic arbitrage. Plugging in the 
BSAi formula for B and its partial derivatives ^, ioTi dsda ^ ^~2"i we can reduce the "PDE" constraint 
into an algebraic restriction on the shape of the implied volatility surface I{t; K, T) 



I\t-K,T) 



fi{t)I{t-K,T)T 



-^ - p{t)uj{t)^/^y/^d2 + '^^-^dld2T 







where t = T — t 

This algebraic restriction is the basis for the specific VGVV models: SRV and LNV, that we describe 
next. 

For SRV we assume square-root implied variance dynamics 

dl^{t) = K{t) [e{t) - /2(t)] dt + 2u;(t)e-''W(^-*)/(t)(iZ(i) 
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If we represent the implied volatility surface in terms of r = T — t and standardized moneyness z{t) = 
'^^ ' 'r- — -, then I(z,t) solves the following quadratic equation 



Z,T 



(1 + K{t)) I\z, r) + (u;2(i)e~2''Wvi-5z) /( 



For LNV we assume log-normal implied variance dynamics 

dl\t) = K{t) [9{t) - l\t)] dt + 2'u;(t)e-''W(^-*)/(t)dZ(i) 

If we represent the implied volatility surface in terms of r = T — t and log relative strike k(t) = In (■^/5(t)), 
then I{k, r) solves the following quadratic equation 



w^it) 



e-^^it)-T^I^{k,T) + fl + K{t)T + u;2(t)e^2ry(t)r^ _ p(^t)y^w{t)e~'^^'^^] l\k,T) 



- v{t) + K9{t)T + 2p(t)y^e"^W"A; + u;2(t)e-2r?{t)r^2 
For both SRV and LNV models we have six time varying stochastic coefficients: 

K{t),e{t),w{t),r,{t),p{t),v{t) 

Given time t values for the six coefficients, the whole implied volatility surface at time t can be found 
as solution to quadratic equations. 

The dynamic calibration procedure treats the six coefficients as a state vector X{t) and it assumes that 
X{t) propagates like a random walk 



where T^x is a diagonal matrix. It also assumes that all implied volatilities are observed with errors 
IID normally distributed with error variance Cg 

y{t) = h{X{t)) + ^ye{t) 

with h{-) denoting the model value (quadratic solution for SRV or LNV) and Sy = Inc^I, with /jv 
denoting an identity matrix of dimension N 

This setup introduces seven auxiliary parameters B that define the covariance matrix of the state and 
the measurement errors. 

When the state propagation and the measurement equation are Gaussian linear, the Kalman filter 
provides efficient forecasts and updates on the mean and covariance of the state and observations. The 
state-propagation equations are Gaussian and linear, but the measurement functions h {X{t)) are not linear 
in the state vector. To handle the non-linearity we employ the unscented Kalman filter. For additional 
details the reader is referred to [■>(')]. 

The procedure was applied successfully on both currency options and equity index options, and compared 
with Heston. 

The comparison with Heston provided the following conclusions: 

• generated half the root mean squared error 
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• explains 4% more variation 

• generated errors with lower serial correlation 

• can be calibrated 100 times faster 

• The whole sample (573 weeks) of implied volatility surfaces can be fitted in about half a second 
(versus about 1 minute for Heston). 

5 Volatility surface based on parametric representations 

Various parametric or semi-parametric representations of the volatility surface have been considered in 
the literature. A recent overview was given in [62]. 

5.1 Polynomial parametrization 

A popular representation was suggested in ['>('>], which proposed that the implied volatility surface is 
modeled as a quadratic function of the moneyness Ai = ln(^/A')/VT 

a{M,T) = bi + b2M + hM^ + b4T + b^MT 

This model was considered for oil markets in [2S], concluding that the model gives only an "average" 
shape, due to its inherent property of assuming the quadratic function of volatility versus moneyness to 
be the same across all maturities. Note that increasing the power of the polynomial volatility function 
(from two to three or higher) does not really offer a solution here, since this volatility function will still 
be the same for all maturities. 

To overcome those problems a semi parametric representation was considered in [28], where they kept 
quadratic parametrization of the volatility function for each maturity T, and approximate the implied 
volatility by a quadratic function which has time dependent coefficients. 

A similar parametrization (but dependent on strike and not moneyness) was considered in [39] under the 
name Practitioner's BlackScholes. It was shown that outperforms some other models in terms of pricing 
error in sample and out of sample. 

Such parametrizations may some certain drawbacks, such as: 

• are not designed to ensure arbitrage-free of the resulting volatility surface 

• the dynamics of the implied volatility surface may not be adequately captured 
We now describe other parametrizations that may be more suitable. 

5.2 Stochastic volatility inspired (SVI) parametrization 

SVI is a practitioner designed parametrization [76, 77]. Very recent papers provide the theoretical frame- 
work and describe its applicability to energy markets [52, 53]. We also note that SVI procedure may be 
employed together with conditions for no vertical and horizontal spread arbitrage, such as in [82]. 

The essence of SVI is that each time slice of the implied volatility surface is fitted separately, such 
that in the logarithmic coordinates the implied variance curve is a hyperbola, and additional constraints 
are imposed that ensure no vertical/ horizontal spread arbitrage opportunities. The hyperbola is chosen 
because it gives the correct asymptotic representation of the variance when log-strike tends to plus or 
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minus infinity: written as a function of ln{^/F), where K is the strike and F is the forward price, and 
time being fixed, the variance tends asymptotically to straight lines when In {^/f) — )• ±00 
The parametrization form is on the implied variance: 

a'^ [x] = v{{m, s,a,b, p} ,x) = a + b I p{x — m) + y {kx — m) + s'' 

where o, 6, p, m, s are parameters which are dependent on the time slice and x = In {^/f). 

We should note that it was recently shown [127] that SVI may not be arbitrage-free in all situations. 
Nevertheless SVI has many advantages such as small computational time, relatively good approximation 
for implied volatilities for strikes deep in- and out-of-the-money. The SVI fit for equity markets is much 
better than for energy markets, for which [53] reported an error of maximum 4-5% for front year and 
respectively 1-2% for long maturities. 

Quasi explicit calibration of SVI is presented in [51], based on dimension reduction for the optimization 
problem. The original calibration procedure is based on matching input market data l(jf'^^'^\ ■_^ yw-, 
which becomes an optimization problem with five variables: a, 6, p, m, s: 



mm > I V 

{a,b,p,m,s} ^ 

1=1 



{m,s,a,b,p} ,ln ( -^ 



E 



2 



The new procedure is based on a change of variables 

X — m 



y 



s 
Focusing on total variance V = vT, the SVI parametrization becomes 









V{y) = 


aT + 5y + /3vV 


+ 1 


re we 


; have used the fohowing 


notations 
















/3 = 


= bsT 












5 -- 


= pbsT 












a = 


= aT 




also 


use the notation Vi = 


^ '^MKT^'^rp 









Therefore, for given m and s, which is transformed into |yj,Vi}, we look for the solution of the 3- 
dimensional problem 



min,^{i/.yj(/^''5'«) (5-1) 

{P,o,a\ 



with the objective functional for reduced dimensionality problem defined by 

N , , s 2 



P{y^,v,}^^^ 5, a) = ^ ?yj a + 5yi + ^^y} + 1 - V- j 



The domain on which to solve the problem is defined as 

^MIN < /3 < 4s 
-/? < (5 < /3 
(4s - /3) < 5 < (4s - /?) 

OiMIN < a < ^MAX 
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For a solution {f3* ,6* ,a*}oi the problem (5.1), we identify the corresponding triplet {a* ,b* , p*} and 
then we solve the 2-dimensional optimization problem 



mm 



{m,s} ^^^ 



{m, s, a*, b*,p*}, In ( -^ 



V 



MKT \ 



Thus the original calibration problem was cast as a combination of distinct 2-parameter optimization 
problem and, respectively, 3-parameter optimization problem. Because the "2+3" procedure is much less 
sensitive to the choice of initial guess, the resulting parameter set is more reliable and stable. For additional 
details the reader is referred to [51]. The SVI parametrization is performed sequentially, expiry by expiry. 
An enhanced procedure was presented in [82] to obtain a satisfactory term structure for SVI, which satisfies 
the no-calendar spread arbitrage in time while preserving the condition of no-strike arbitrage. 

5.3 Entropy based parametrization 

Entropic calibrations have been considered by a number of authors. It was done for risk-neutral terminal 
price distribution, implied volatility function and the option pricing function. 

An algorithm that yields an arbitrage-free diffusion process by minimizing the relative entropy distance 
to a prior diffusion is described in [8]. This results in an efficient calibration method that can match 
an arbitrary number of option prices to any desired degree of accuracy. The algorithm can be used to 
interpolate, both in strike and expiration date, between implied volatilities of traded options. 

Entropy maximization is employed in [31] and [30] to construct the implied risk- neutral probability 
density function for the terminal price of the underlying asset. The advantage of such an entropic pricing 
method is that it does not rely on the use of superfluous parameters, and thus avoids the issue of over 
fitting altogether. Furthermore, the methodology is flexible and universal in the sense that it can be 
applied to a wider range of calibration situations. 

Most of the entropy-based calibration methodologies adopted in financial modeling, whether they are 
used for relative entropy minimization or for absolute entropy maximization, rely on the use of the logarith- 
mic entropy measure of Shannon and Wiener. One drawback in the use of logarithmic entropy measures is 
that if the only source of information used to maximize entropy is the market prices of the vanilla options, 
then the resulting density function is necessarily of exponential form. On the other hand, empirical studies 
indicate that the tail distributions of asset prices obey power laws of the Zipf-Mandelbrot type [30] . Thus 
we would like to employ entropies that may recover power law distribution, such as Renyi entropy [30] 

Maximization of Renyi entropy is employed to obtain arbitrage-free interpolation. The underlying 
theoretical idea is that, irrespective of the nature of the underlying price process, the gamma associated 
with a European-style vanilla option always defines a probability density function of the spot price implied 
by the existence of the prices for option contracts. There is a one-to-one correspondence between the 
pricing formula for vanilla options and the associated gamma. Therefore, given option gamma we can 
unambiguously recover the corresponding option pricing formula. 

We present here an overview of the approach presented in [30] 

Given strikes Kj, j = 1...M, corresponding for input market prices, maximizing the Renyi entropy yields 
a density function of the form: 







p{x) = I A + /3ox + V /3, (x - Kj)^ I (5.2) 
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The parameters a, A, Pq, ..., (3^ are calibrated by matching the input prices to the prices computed using 
the density function (5.2). 

We exemphfy for the case of cah options. For each j = 1...M, we have to impose the matching condition 
to market price C^ 

m=l ^ ^ 

where 




We also impose the normalization condition 

7 _ M 

J p{x)dx = 1 =^ ^ Yl ^i [^i (^)] ^ 



iX=Ki 



j=0 



Since the density function is explicitly given, is straightforward to use for calibration additional option 
types, such as digitals or variance swaps. 

The result is described in [30] as leading to the power-law distributions often observed in the market. 
By construction, the input data are calibrated with a minimum number of parameters, in an efficient 
manner. The procedure allows for accurate recovery of tail distribution of the underlying asset implied by 
the prices of the derivatives. One disadvantage is that the input values are supposed to be arbitrage free, 
otherwise the algorithm will fail. It is possible to enhance the algorithm to handle inputs with arbitrage, 
but the resulting algorithm will lose some of the highly efficient characteristics, since now we need to solve 
systems of equations in a least square sense 

5.4 Parametrization using weighted shifted lognormal distributions 

A weighted sum of interpolation functions taken in a parametric family is considered in the practitioner 
papers [26, 27] to generate a surface without arbitrage in time and in space, while remaining as closely as 
possible to market data. Each function in the family is required to satisfy the no-free lunch constraints, 
specified later, in such way that they are preserved in the weighted sum. 

In this parametric model, the price of a vanilla option price of strike K and maturity T is estimated at 
time to = by the weighted sum of N functionals 

N 



Yai{T)Fi{toSo,P{Ty,K,T) 



i=l 



with ai{T) weights and P{t) = B (0,t) the zero coupon bond price. 

Several families Fi can satisfy the No-Free-Lunch constraints, for instance a sum of lognormal distri- 
butions, but in order to match a wide variety of volatility surfaces the model has to produce prices that 
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lead to risk-neutral pdf of the asset prices with a pronounced skew. If all the densities are centered in 
the log-space around the forward value, one recovers the no-arbitrage forward pricing condition but the 
resulting pricing density will not display skew. However, centering the different normal densities around 
different locations (found appropriately) and constraining the weights to be positive, we can recover the 
skew. Since we can always convert a density into call prices, we can then convert a mixture of normal 
densities into a linear combination of BSM formulae. 

Therefore, we can achieve that goal with a sum of shifted log-normal distributions, that is, using the 
BSAi formula with shifted strike (modified by the parameters) as an interpolation function 

Fi (to, So, P{T); K, T) = CalhsM (to, So, P{T), K (1 + ^x,{T)) , T, a. 



with i^denoting adjusted strike. 

We note that the value of strike is adjusted only if we apply the procedure for equity markets, in which 
case it becomes 

k{K,t) =K + D{0,t) 

with D{0,t) is the compounded sum of discrete dividends between and t. 

The no-arbitrage theory imposes time and space constraints on market prices. Hence, the time dependent 
parameters ai{t) and Hi{t) are used to recover the time structure of the volatility surface. It is argued that 
it sufficient to use a parsimonious representation of the form 

N 



ai{t) 
fit,|3^) 



V- «M «° 



2 



1+(1 + ^ 



2 



Making the weights and the shift parameter time-dependent to fit a large class of volatility surfaces 
leads to the following no- free lunch constraints, for any time t 

• flj > to get convexity of the price function 

• X^j=i o,i{t) = 1 to get a normalized risk- neutral probability 

• X^i=i o,i{t)fii{t) = 1 to keep the martingale property of the induced risk-neutral pdf 

• /ii(t) to get non-degenerate functions 

The model being invariant when multiplying all the terms a^ with the same factor, we impose the nor- 
malization constraint 

N 

to avoid the possibility of obtaining different parameter sets which nevertheless yield the same model. 

Given the A^ parameters and assuming a constant volatility ai{t) = cr,-', there are AN — 2 free parameters 
for theA^-function model since we can use the constraints to express a^ and, respectively, a^fii in terms of 
{oi'}_i ly ^^'^ {^^^} -1 N (^^^ ^^^° Appendix A of [27]). 
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As such, this model does not ahow for the control of the long term volatility surface. Therefore, for the 
model to be complete we specify the time-dependent volatility such that it captures the term structure of 
the implied volatility surface: 

ai{t)=-fie-'^' + dif{t,bi) 

Thus we need to solve a IN — 2 optimization problem. This is done in [26, 27] using a global optimizer 
of Differential Evolution type. 

6 Volatility surface based on nonparametric representations, smoothing 
and interpolation 

This broad set of procedures may be divided into several categories. 

6.1 Arbitrage-free algorithms 

Interpolation techniques to recover a globally arbitrage-free call price function have been suggested in 
various papers, e.g., [95, 143]. A crucial requirement for these algorithms to work that the data to 
be interpolated are arbitrage-free from the beginning. [95] proposes an interpolation procedure based on 
piecewise convex polynomials mimicking the BSM pricing formula. The resulting estimate of the call price 
function is globally arbitrage-free and so is the volatility smile computed from it. In a second step, the 
total (implied) variance is interpolated linearly along strikes. Cubic B-spline interpolation was employed 
by [14.]], with interpolation performed on option prices, and the shape restrictions in interpolated curves 
was imposed by the use of semi-smooth equations minimizing the distance between the implied risk neutral 
density and a prior approximation. 

Instead of smoothing prices, [20] suggests to directly smooth implied volatility parametrization by means 
of constrained local quadratic polynomials. Let us consider that we have M expiries {Tj} and N strikes 
{xi} , while the market data is denoted by {cf (Tj)^ 

Two approaches are considered: 

• each maturity is treated separately 

• all maturities are included in the cost functional to minimize 

The first case implies minimization of the following (local) least squares criterion at each expiry Tj , j = l...NT 



N 

min y |af ^^(T,) - a^^^ - a? {xi - x) - 4'^ {xi - xf] y} 






h 



where /C is a kernel function, typically a symmetric density function with compact support. 
One example is the Epanechnikov kernel 

/C(n) = 0.75(l-n2) 1 [|n| < 1] 

with 1{A) denoting the indicator function for a set A and h is the bandwidth which governs the trade-off 
between bias and variance. 

The optimization problem for the second approach is 



A'l N 









hx 



hx 



hx 
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with defined as 



-ai~^ {Tj -T)- af {xi - xf - af^ {xi - x) {Tj - T) 

The approach yields a volatihty surface that respects the convexity conditions, but neglects the condi- 
tions on call spreads and the general price bounds. Therefore the surface may not be fully arbitrage-free. 
However, since convexity violations and calendar arbitrage are by far the most virulent instances of arbi- 
trage in observed implied volatility data, the surfaces will be acceptable in most cases. 

The approach in [02] is based on cubic spline smoothing of option prices rather than on interpolation. 
Therefore, the input data does not have to be arbitrage-free. It employs cubic splines, with constraints 
specifically added to the minimization problem in order to ensure that there is no arbitrage. A potential 
drawback for this approach is the fact that the call price function is approximated by cubic polynomials. 
This can turn out to be disadvantageous, since the pricing function is not in the domain of polynomials 
functions. It is remedied by the choice of a sufficiently dense grid in the strike dimension. 

Instead of cubic splines, [102] employs constrained smoothing B-splines. This approach permits to im- 
pose monotonicity and convexity in the smoothed curve, and also through additional pointwise constraints. 
According to the author, the methodology has some apparent advantages on competing methodologies. It 
allows to impose directly the shape restrictions of no-arbitrage in the format of the curve, and is robust the 
aberrant observations. Robustness to outliers is tested by comparing the methodology against smoothing 
spline. Local Polynomial Smoothing and Nadaraya- Watson Regression. The result shows that Smoothing 
Spline generates an increasing and non-convex curve, while the Nadaraya- Watson and Local Polynomial 
approaches are affected by the more extreme points, generating slightly non convex curves. 

It is mentioned in [117] that a large drawback of bi-cubic spline or B-spline models is that they require the 
knots to be placed on a rectangular grid. Correspondingly, it considers instead a thin-spline representation, 
allowing arbitrarily placed knots. This leads to a more complex representation at shorter maturities while 
preventing overfitting. 

Thin-spline representation of implied volatility surface was also considered in [29] and section 2.4 of 
[96], where it was used to obtain a pre-smoothed surface that will be eventually used as starting point for 
building a local volatility surface. 

An efficient procedure was shown in [109] for constructing the volatility surface using generic volatility 
parametrization for each expiry, with no-arbitrage conditions in space and time being added as constraints, 
while a regularization term was added to the calibrating functional based on the difference between market 
implied volatilities and, respectively, volatilities given by parametrization. Bid-ask spread is also included 
in the setup. The resulting optimization problem has a lot of sparsity/structure, characteristics that were 
exploited for obtaining a good fit in less than a second 

6.2 Remarks on spline interpolation 

The following splines are usually employed to interpolate implied volatilities 

• Regular cubic splines 

• Cubic B-splines 

• Thin splines 

Certain criteria (such as arbitrage free etc) have to be met, and relevant papers were described in the 
previous section . Here we just refer to several generic articles on spline interpolation. 
[145] describes an approach that yields monotonic cubic spline interpolation. 
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Although somewhat more comphcated to implement, B-splines may be preferred to cubic splines, due 
to its robustness to bad data, and ability to preserve monotonicity and convexity. A recent paper [99] 
describes a computationally efficient approach for cubic B-splines. 

A possible alternative is the thin-plate spline, which gets its name from the physical process of bending 
a thin plate of metal. A thin plate spline is the natural two-dimensional generalization of the cubic spline, 
in that it is the surface of minimal curvature through a given set of two-dimensional data points. 

6.3 Remarks on interpolation in time 

In some situations we need to perform interpolation in time. While at a first glance it may seem straightfor- 
ward, special care has to be employed to ensure that the result still satisfies practical arbitrage conditions. 
For example, one should expect that there is no calendar spread arbitrage [34, 52, 76, 124] 

One common approach is to perform linear interpolation in variance. A variant of it, denoted "total 
variance interpolation", is described in [37]. 

6.4 Interpolation based on fully implicit finite difference discretization of Dupire forward 
PDE 

We present an approach described in [6, 7, 91], based on fully implicit finite difference discretization of 
Dupire forward PDE. 

We start from the Dupire forward PDE in time-strike space 



dc 1 r / , m2 d c 



Let us consider that we have the following time grid = Iq < ti < ... < tj\[ and define Atj = tj+i — tj 
A discrete (in time) version of the forward equation is 

c{ti+i,k)-c{ti,k) 1 2 d^c 
A^^ = 2[^(*-^)] 9p('^+^'^) 

This is similar to an implicit finite difference step. It can be rewritten as 



l-^At,[aiU,k)f^ 



c{ti+i,k) = c{ti,k) (6.1) 



Let us consider that the volatility function is piecewise defined on the time interval ti < t < ij+i and 
we denote by i'i{k) the corresponding functions 

i^i{k) = a{t, k) for U < t < tj+i 

Using (6.1) we can construct European (call) option prices for all discrete time points for a given a set 
of volatility functions {i^iik)}^^^ ^ by recursively solving the forward system 



l-^AU[aiU,k)f^ 



c{U+i,k) = c{ti,k) (6.2) 

c{0,k) = [S{0)-k] + 
Let us discretize the strike space as Km in = kQ < ki < ... < kM = kMAX 



23 



By replacing the differential operator 9'^/dk'^ by the central difference operator 

2 „,, . 2 



Skkf{k) 



{kj - kj_i) {kj+i - /cj_i) ^ {kj - kj-i) {kj+i - kj) ^ 



+ - 



:m 



{k,+,-k,){k,+l-k,.ly'^^' 

we get the following finite difference scheme system 



1 2 

1 - -Atj [ui {k)] 6kk 



c{ti+i,k) 
c(0,/c) 



c{ti,k) 
[5(0) - k]- 



(6.3) 



The matrix of the system (6.3) is tridiagonal and shown in [6] to be diagonally dominant, which allows 
for a well behaved matrix that can be solved efficiently using Thomas algorithm [132]. Thus we can directly 
obtain the European option prices if we know the expressions for {fj(A;)}^^^ ^. 

This suggests that we can use a bootstrapping procedure, considering that the volatility functions are 
defined as piecewise constant 

Let us first introduce the notations for market data. We consider that we have a set of discrete option 

is the set of strikes for expiry 

and, 



quotes \^c^^^'^ [ti,Ki^p)\, where {tj} are the expiries and {Ki^p} ^-^ nkH) 
ti. 



We should note that we may have different strikes for different expiries, and that {Ki^p} -^ NK(i) 
respectively, {kj} represent different quantities 

Then the piecewise constant volatility functions are denoted as 



u,{k) 



ai^p for Ki,p <k< Ki., 



P+i 



Thus the algorithm consists of solving an optimization problem at each expiry time, namely 
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{o-i,l,---,0-i,NK(i)} 



'NK{i) 

Y, {c{U,K,,p)-c^'''''{U,K,,p)y 

p=l 



(6.4) 



We remark that, when solving (6.4) by some optimization procedure, one needs to solve only one 
tridiagonal matrix system for each optimization iteration. 

Regarding interpolation in time, two approaches are proposed in [6]. The first one is based on the 
formula 

I,. ... ..,2 d^ 



i--{t-u) w, {k)Y ^ 



c{ti+i,k) = c{ti,k) for ti <t <ti+i 
while the second one is a generalization of (6.5) 



1 - i (r(t) - u) [n {k)f ^ 



c{ti+i,k) = c{ti,k) for ti <t <ti+i 



(6.5) 



(6.6) 



where T{t) is a function that satisfies the conditions T{ti) = ti and T'{t) < 
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It is shown in [G] that option prices generated by (6.2)and (6.5) and, respectively, by (6.2) and (6.6)are 
consistent with the absence of arbitrage in the sense that , for any pair (t, k) we have 

dc 

7 Adjusting inputs to avoid arbitrage 

Various papers have tackled the problem of finding conditions that may be necessary/and or sufficient 
to ensure that prices/vols are free of arbitrage [34, 46, 49, 89, 113, 125]. If one wants to adjust the set 
of input prices/vols to avoid arbitrage, several approaches have been described in the literature. For 
example, [3] presents a relatively simple method to adjust implied volatilities, such that the resulting set 
is both arbitrage free and also closest to the initial values (in a least-squares sense). Another algorithm 
is presented in section 8.3 of [32]. We present in detail the algorithm from [34], based on the observation 
that the absence of call spread, butterfly spread and calendar spread arbitrages is sufficient to exclude all 
static arbitrages from a set of option price quotes across strikes and maturities. 

7.1 Carr and Madan algorithm 

The main idea is as follows: given input market prices and corresponding bid ask spreads, we start from 
the price corresponding to first expiry ATM and adjust the prices for that expiry. We continue to the 
next expiry and we make sure that arbitrage constraints are satisfied both in time and strike space, while 
adjusting within the bid ask spread. 

We present first the arbitrage constraints from [34], using notations from there. Let Cij denote the given 
quote for a call of strike Ki and maturity Tj. We suppose that the N strikes {Ki} form an increasing and 
positive sequence as do the M maturities {Tj}. Without any loss of generality, we suppose that interest 
rates and dividends are zero over the period ending at the longest maturity. 

We augment the provided call quotes with quotes for calls of strike Kq = 0. For each maturity, these 
additional quotes are taken to be equal to the current spot price Sq. We also take the prices at maturity 
Tq = to be (5o — Ki) . This gives us the augmented matrix of prices Cij, with indices i = O..A^ and 
j = 1...M . 

For each j > we define the following quantities: 

"^''^ - K, - K,., 

Qoj = 

For each i > 0, Qij is the cost of a vertical spread which by definition is long y{Ki-Ki-i) calls of strike 
Ki-i and short Y(A'i-A'i_i) calls of strike Ki. A graph of the payoff from this position against the terminal 
stock price indicates that this payoff is bounded below by zero and above by one. 

We therefore require for our first test that 

< Qij <1, i = 1...N, j = 1...M (7.1) 

Next, for each j > 0, we define the following quantities: 

oc A^ Ki+i-Ki_i Ki-Ki_i 
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For each i > BSprij is the cost of a butterfly spread which by definition is long the call struck at 
Ki-i, short {^i+i-^i-i)/(Ki+i-Ki) calls struck at Ki, and long {^i-^i-i)/{KiJ^i-Ki) calls struck at -ffj+i. A 
graph indicates that the butterfly spread payoff is non-negative and hence our second test requires that 



Ki+i - Ki_i ^ . Ki- Ki_i 
Equivalently, we require that 






We define 



Ki - Aj_i Aj+i - Ki 

We may interpret each qij as the marginal risk-neutral probability that the stock price at maturity Tj 
equals Ki. 

For future use, we associate with each maturity a risk-neutral probability measure defined by 



liK) = Yl 1^^^ 



A third test on the provided call quotes requires that for each discrete strike Ki,i > 0, and each discrete 
maturity Tj,j > 0,we have 

Qj+i - Cj > (7.3) 

The left-hand side of (7.3) is the cost of a calendar spread consisting of long one call of maturity Tj^i 
and short one call of maturity Tj, with both calls struck at Ki. Hence, our third test requires that calendar 
spreads comprised of adjacent maturity calls are not negatively priced at each maturity. 

We now conclude, following [34] , the discussion on the 3 arbitrage constraints (7.1)(7.2)(7.3). 

As the call pricing functions are linear interpolations of the provided quotes, we have that at each 
maturity Tj, calendar spreads are not negatively priced for the continuum of strikes K > 0. Since all 
convex payoffs may be represented as portfolios of calls with non-negative weights, it follows that all 
convex functions (j){S) are priced higher when promised at Tj+i than when they are promised at Tj. In 
turn, this ordering implies that the risk-neutral probability measures Q constructed above are increasing 
in the convex order with respect to the index j. This implies that there exists a martingale measure which 
is consistent with the call quotes and which is defined on some filtration that includes at least the stock 
price and time. Finally, it follows that the provided call quotes are free of static arbitrage by standard 
results in arbitrage pricing theory. 

8 Characteristics of volatility surface 

Many recent papers have studied various characteristics of volatility surface: 

• the static and dynamic properties of the implied volatility surface must exhibit within an arbitrage- 
free framework 

• implied volatility calculations in a (local) stochastic volatility environment, which may also include 
jumps or even Levy processes. 
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• the behavior of impUed volatihty in Hmiting cases, such as extreme strikes, short and large maturities, 
etc. 

For completion we include a list of relevant papers: [11, 14, 21, 38, 122, 46, 48, 49, 50, 57, 63, 70, 69, 66, 
65, 64, 67, 68, 71, 73, 75, 77, 78, 81, 86, 87, 93, 101, 103, 104, 106, 112, 111, 115, 118, 119, 123, 125, 127, 
128, 133, 134, 135, 137, 140, 141, 16, 17, 19, 22, 23, 24, 25, 10, 58, 72, 74, 92, 97, 136, 144, 15] 

The constructed volatility surface may also need to take into account the expected behavior of the 
volatility surface outside the core region. The core region is defined as the region of strikes for equity 
markets, moneyness levels for Commodity markets, deltas for FX markets for which we have observable 
market data. From a theoretical point of view, this behavior may be described by the asymptotics of 
implied volatility, while from a practical point of view this corresponds to smile extrapolation. 

8.1 Asymptotics of implied volatility 

Concerning the dependence with respect to strike, some major theoretical results are known in a model- 
independent framework. [103] related the extreme strike slopes of the implied volatility to the critical 
moments of the underlying through the moment formula: let a{t, x) denote the implied volatility of a 
European Call option with maturity t and strike K = Fqc^, then 

limsup^^^^^^=VK(t)-l) (8.1) 

a;— > oo X 

where ip{u) = 2 — 4 ( \/v?+u — u I and u* (t) = sup {u> l;E [F'^{t)]} is the critical moment of the 

underlying price F = (F(t))j>Q. An analogous formula holds for the left part of the smile, namely when 
X —7- — oo. This result was sharpened in [14, 15], relating the left hand side of (8.1) to the tail asymptotics 
of the distribution of F(t). 

In the stochastic volatility framework this formula was applied by [">] and [97], to mention but a few. 

The study of short- (resp. long-) time asymptotics of the implied volatility is motivated by the research 
of efficient calibration strategies to the market smile at short (resp. long) maturities. Short time results 
have been obtained in [111, 66, 65, 64, 21], while some other works provide insights on the large-time 
behavior, as done by [141] in a general setting, [97] for affine stochastic volatility models or [07] for Heston 
model. 

8.2 Smile extrapolation 

It is argued in the practitioner paper [J 3] that a successful smile extrapolation method should deliver 
arbitrage-free prices for the vanilla options, i.e., the option prices must be convex functions of strike, 
and stay within certain bounds. In addition, the extrapolation method should ideally have the following 
properties: 

1. It should reprice all observed vanilla options correctly. 

2. The PDF, CDF and vanilla option prices should be easy to compute. 

3. The method should not generate unrealistically fat tails, and if possible, it should allow us to control 
how fat the tails are. 

4. It should be robust and flexible enough to use with a wide variety of different implied volatility 
surfaces. 
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5. It should be easy and fast to initialize for a given smile. 

The paper describes two commonly used methods which do not satisfy the above wish list. The first one 
is to use some interpolation within the region of observed prices, and just set the implied volatility to be a 
constant outside of this region. This method is flawed as it introduces unstable behavior at the boundary 
between the smile and the flat volatility, yielding unrealistically narrow tails at extreme strikes. 

The second approach is to use the same parametric form for the implied volatility derived from a model, 
e.g. SABR, both inside and outside the core region. There are several problems with this method. It 
gives us little control over the distribution; indeed this approach often leads to excessively fat tails, which 
can lead to risk neutral distributions that have unrealistically large probabilities of extreme movements, 
and have moment explosions that lead to inflnite prices, even for simple products. If the methodology is 
dependent on usage of an asymptotic expansion, the expansion may become less accurate at strikes away 
from the money, leading to concave option prices, or equivalently negative PDFs, even at modestly low 
strikes. Furthermore, there is no guarantee that this functional form will lead to arbitrage free prices for 
very large and small strikes. 

That is why [13] propose to separate the interpolation and extrapolation methods. 

Their method works as follows. A core region of observability, inside which we may use any standard 
smile interpolation method, is deflned first: i^_ < K < K^. Outside of this region the extrapolation is 
done by employing a simple analytical formula for the option prices, that has the following characteristics: 



• 



For very low strikes region, the formula-based put prices will go towards zero as the strike goes to 
zero, while remaining convex. 



• 



For very high strikes region, the formula-based call prices will go towards zero as strike goes to 
infinity, while remaining convex. 

Each of these formulas is parametrized so that we can match the option price as well as its first two 
derivatives at the corresponding boundary with the core region. The methodology is also able to retain a 
measure of control over the form of the tails at extreme strikes. 

The following functional form for the extrapolation of put and, respectively, call prices was described 
as parsimonious yet effective: 

Put{K) = Ki" exp [ai + hK + ciK^] 



Call{K) =K-'' exp 



b2 C2 

do -\ TT 

K K^ 



We fix fi>l, which ensures that the price is zero at zero strike, and there is no probability for the 
underlying to be zero at maturity. Alternatively, we can choose fi to reflect our view of the fatness of the 
tail of the risk neutral distribution. It is easy to check that this extrapolation generates a distribution 
where the ?TT,-th negative moment is finite for m < 1 — n and infinite for m > 1 — /_i. 

We fix i^ > to ensure that the call price approaches zero at large enough strikes. Our choice of controls 
the fatness of the tail; the m-th moment will be finite if tti < z^ — 1 and infinite if m > v — 1. 

The condition for matching the price and its first two derivatives at K^ and, respectively, at K^ yields 
a set of linear equations for the parameters ai,6i,ci and, respectively, for a2,b2,C2 

9 Remarks on numerical calibration 

The calibration procedure consists of finding the set of parameters (defining the volatility surface) that 
minimize the error between model output and market data, while satisfying some additional constraints 
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such as "no arbitrage in strike and time", smoothness, etc. This chapter provides some details regarding 
the practical aspects of numerical calibration. 

Let us start by making some notations: we consider that we have Mexpiries JT"'} and that for each 
maturity T^^' we have N[j] calibrating instruments, with strikes Kij, for which market data is given (as 
prices or implied volatilities). The bid and ask values are denoted by Bid [i,T^^'\ and Ask [i^T^^n 

9.1 Calibration function 

The calibration function is defined in different ways 

If we perform "all-at-once" calibration, then the calibration function is constructed as 

M N[j] 

^ = ^Y^ Wi,j ModelOutput (i, T'^^A - MarketData (i, T^^A 
j=i i=i 

where Wij are weights and ||-|| denotes a generic norm 

If we perform sequential calibration, one expiry at the time, than the calibration functional for each 
expiry will be given as 

^ [j] A Y^ Wij ModelOutput (i, T^^A - MarketData (i, T^J^ ' 

i=l 

If we use a local optimizer, then we might need to add a regularization term. The regularization term 
the most commonly considered in the literature is of Tikhonov type, e.g., [116, 44, 110, 2]. However, since 
this feature is primarily employed to ensure that the minimizer does not get stuck in a local minimum, 
this additional term is usually not needed if we use either a global optimizer or a hybrid (combination of 
global and local) optimizer. 

9.2 Constructing the weights in the calibration functional 

The weights Wij can be selected following various procedures detailed in [26, 27, 41, 47, 105] chapter 13 
of [42], to mention but a few. 

Practitioners usually compute the weights (see [J 7]) as inverse proportional to 

• the square of the bid-ask spreads, to give more importance to the most liquid options. 

• the square of the BSM Vegas (roughly equivalent to using implied volatilities, as explained below). 

[105] asserts that it is statistically optimal (minimal variance of the estimates) to choose the weights as 
the inverse of the variance of the residuals, which is then considered to be proportional to the inverse of 
squared bid-ask spread. 

Another practitioner paper [26] considers a combination of the 2 approaches and this is our preferred 
methodology. 

It is known that at least for the options that are not too far from the money, the bid-ask spreads is of 
order of tens of basis points. This suggests that it might be better to minimize the differences of implied 
volatilities instead of those of the option prices, in order to have errors proportional to bid-ask spreads 
and to have better scaling of the cost functional. However, this method involves additional computational 
cost. A reasonable approximation is to minimize the square differences of option prices weighted by the 
BSM Vegas evaluated at the implied volatilities of the market option prices. 
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The starting point is given by setting the weights as 

1 



w^ 



«j 



|Bi(i(i,rO)) - Ask{i,Tiiy 



To simphfy the notations for the remainder of the chapter we denote the model price by ModP and the 
market price by MktP 

We can approximate the difference in prices as foUows: 

ModP (z,rO)) - MktP (^,^0•)) ^ ^[^''^/J^^^^'^)] {afiO^ {^,T(^)) - a%^^ (z,rO) 

where crfy {i,T^^>^ and crfy {i,T^^'^ are the imphed vols corresponding to model and, respectively, 
market prices for strikes Ki^j and maturities T^^'. 

We should note that this series approximation may be continued to a higher order if necessary, with a 
very small additional computational cost. 

Using the following expression for BSM Vegas evaluated at the implied volatility crfy^'^ (i, T"') of the 
market option prices 

we obtain 

Thus we can switch from difference of implied volatilities to difference of prices. For example, for all- 
at-once calibration that is based on minimization of root mean squared error (RMSE) , the corresponding 
calibration functional can be written as 

M N[j] 

^ = ^ ^ tiJij ModelOutput (i, t'-^A - MarketData (i, T^^^ 

where the weights Wij are defined as 

1 f 1 ^ 

Wi 



^«j 



Bid{i,TU))-Ask {i,TU))\ \Vega{af(/^^ (i,rO))) 



To avoid overweighting of options very far from the money we need introduce an upper limit for the 
weights. 

9.3 Selection of numerical optimization procedure 

It is quite likely that the calibration function for the volatility surface may exhibit several local (and 
perhaps global) minima, making standard optimization techniques somewhat unqualified according to 
[40]. Gradient based optimizers, for example, are likely to get stuck in a local minimum which may also be 
strongly dependent on the initial parameter guess. While this situation (multiple local minima) may be 
less common for equity markets , it is our experience that such characteristics are quite common for FX 
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and Commodities markets. Thus global/hybrid optimization algorithms are our preferred optimization 
methods in conjunction with volatility surface construction. 

Our favorite global optimization algorithm is based on Differential Evolution [121]. In various flavors it 
was shown to outperform all other global optimization algorithms when solving benchmark problems (un- 
constrained, bounded or constrained optimization). Various papers and presentations described successful 
calibrations done with Differential Evolution in finance: [2G, 27, 54, 142, 40], to mention but a few. 

Here is a short description of the procedure. Consider a search space Q and a continuous function 
G : — 7- M to be minimized on Q. An evolutionary algorithm with objective function G is based on 
the evolution of a population of candidate solutions, denoted by X^ = {0Jj}-_-, j^- The basic idea is to 
"evolve" the population through cycles of modification (mutation) and selection in order to improve the 
performance of its individuals. At each iteration the population undergoes three transformations: 

During the mutation stage, individuals undergo independent random transformations, as if performing 
independent random walks in 0, resulting in a randomly modified population V^ . In the crossover 
stage, pairs of individuals are chosen from the population to "reproduce": each pair gives birth to a new 
individual, which is then added to the population. This new population, denoted W^ , is now evaluated 
using the objective function G{-). Elements of the population are now selected for survival according to 
their fitness: those with a lower value of G have a higher probability of being selected. The A^ individuals 
thus selected then form the new population X^^-^. The role of mutation is to explore the parameter space 
and the optimization is done through selection. 

On the downside, global optimization techniques are generally more time consuming than gradient 
based local optimizers. Therefore, we employ a hybrid optimization procedure of 2 stages which combines 
the strengths of both approaches. First we run a global optimizer such as Differential Evolution for a 
small number of iterations, to arrive in a "neighborhood" of a global minimum. In the second stage we 
run a gradient-based local optimizer, using as initial guess the output from the global optimizer, which 
should converge much quite fast since the initial guess is assumed to be in the "correct neighborhood". An 
excellent resource for selecting a suitable local optimizer can be found at [ I ] 

We should also mention that a very impressive computational speedup (as well as reducing number of 
necessary optimization iterations) can be achieved if the gradient of the cost functional is computed using 
Adjoint method in conjunction with Automatic, or Algorithmic, Differentiation (usually termed AD). Let 
us exemplify the computational savings. Let us assume that the calibration functional depends on P 
parameters, and that the computational cost for computing one instance of the calibration functional is T 
time units. The combination between adjoint and AD methodology is theoretically guaranteed to produce 
the gradient of of the calibration functional (namely all P sensitivities with respect to parameters) in 
a computational time that is not larger than 4-5 times the original time T for computing one instance 
of the calibration functional. The gradient is also very accurate, up to machine precision, and thus we 
eliminate any approximation error that may come from using finite difference to compute the gradient. 
The local optimizer can then run much more efficiently if the gradient of the calibration functional is 
provided explicitly. For additional details on adjoint plus AD the reader is referred to [90] 

10 Conclusion 

We have surveyed various methodologies for constructing the implied volatility surface. We have also 
discussed related topics which may contribute to the successful construction of volatility surface in practice: 
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conditions for arbitrage/non arbitrage in both strike and time, how to perform extrapolation outside the 
core region, choice of cahbrating functional and selection of optimization algorithms. 
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